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ABSTRACT 

An implicit, finite volume code for solving two-dimensional, compressible turbu- 
lent flows is described. Second-order upwind differencing of the inviscid terms of 
the equations is used to enhance stability and accuracy. A diagonal form of the 
implicit algorithm is used to improve efficiency. Several zero- and two-equation 
turbulence models are incorporated to study their impact on overall flow-modeling 
accuracy. Applications to several external and internal flows are discussed. 

1 . INTRODUCTION 

A computer code for solving the time-dependent, Reynolds-averaged, compressible 
Navier— Stokes equations is described. The code solves the two-dimensional (2-D) planar 
and axisymmetric equations in arbitrary, curvilinear coordinates utilizing the finite 
volume approach. The numerical algorithm (refs. 1 and 2) is an approximately factored 
ADI method based on the work of Beam and Warming (ref. 3) and Briley and McDonald 
(ref. 4). It utilizes second-order upwind differencing techniques for the inviscid 
terms of the equations and a diagonal form of the implicit operators. These proce- 
dures are designed to enhance the overall accuracy, stability, and efficiency of the 
algorithm. Rapid convergence to steady-state solutions is achieved through the use 
of spatially varying time steps. A variety of inviscid and viscous boundary condi- 
tions are incorporated. These include subsonic and supersonic inflow and outflow 
conditions, which make use of the method of characteristics, and inviscid and viscous 
solid-wall boundary conditions. 

Various turbulence models are integrated into the code to compare, evaluate, and 
ultimately improve model performance (ref. 5). The models include a zero-equation 
(Cebeci-Smith) model and three two-equation models. The models use the procedure of 
integration to the wall where no-slip conditions at a solid surface are used. 

Although in an early stage of development, the option of using wall functions (or 
law-of-the-wall boundary conditions) in place of no-slip conditions is also included. 

In addition to the above features, the code uses a self-contained, algebraic, 
grid-generation routine, a special data-loading system designed to facilitate implicit 


1 



integration across interior mesh boundaries, a restart capability, and several types 
of input, output, and plotting routines. 

In the following discussion, we open with a description of the theoretical basis 
of the numerical methods and turbulence models. We then discuss several features of 
the code which includes a summary of subroutines. Finally, sample results of invis- 
cid and viscous calculations are presented. 


2 . THEORETICAL BACKGROUND 


2.1 Governing Equations 

The basic equations governing the flow of a viscous compressible fluid are the 
Navier— Stokes equations. These equations, expressed in terms of a stationary two- 
dimensional, curvilinear coordinate system and utilizing vector-dyadic notation, are 
written in strong conservation form as follows: 
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In these equations, £ and ri are the computational coordinates, t is the time, 

^ ~ 3/3t, 3 = 3/3a are partial derivatives, and the summation convention is used 

t ot 

on the repeated subscript a — fjjp. The computational coordinates are related to 
Cartesian coordinates through the coordinate transformation 


x = x(£,ti) , y = yU.n) 


( 2 ) 


The geometric variables in equation (1) are the differential volume element u 
(actually an area in 2-D) and the surface element vectors (actually line elements 

in 2-D), which are defined below. 


u - Vn - Vs 

^5 = ^ s 5x + ^ s Sy = iy n “ Jx n ' 

- ts nx + T s ny “ ^ + 


2 



where 3! and j" are unit Cartesian base vectors and x^ *= S^x, x p = 9 p x, etc * ( see 
fig. 1). An important property of the surface vectors is the metric identity 

3 A - 3 A + 3 A * ° (4) 

which permits ? to be taken either inside or outside the derivatives in operations 
involving the gradient operator i.e., 

^ ° f = ^ ° 3 a f = i 3 a S a o f (5) 

where the symbol o represents any vector product operator. 

The physical variables appearing in equation (1) are the density p, the fluid 
velocity vector u = + jv, the total specific energy E - e -f u • u/2, the specific 

internal energy e, and the pressure p. The variables U and 8F are the fluid-state 
vector and flux tensor (or dyadic), respectively. The flux tensor is split into 

, x , and q appearing in 

are, respectively, the unit tensor, the viscous stress tensor, and the heat flux 
vector, which may be written (using Stokes hypothesis) as 

t" = -y v (^u + uV - • u) 

= ‘ ~ + Wh “ 3 3 3 ^] 

q = -p e ve - - — S 6 3 6 e 

where is the conjugate of Vit For laminar flows, y v = u and y e = yy/Pr are 

the molecular viscosity and conductivity coefficients, where Pr = C p y/K is the 
Prandtl number. It is assumed that the fluid is a perfect gas governed by the equa- 
tion of state 

p = (y - l)pe = pc 2 /y , e = C y T , y = c p / c v < 7 ) 

where y is the ratio of specific heats, c is the sound speed, and T is 
temperature. 

For turbulent flows, the basic equations are time- or ensemble -averaged using 
mass-weighted averaging. The resulting equations are called the Reynolds-averaged 
Navier-Stokes equations and can be cast in a form essentially identical with equa- 
tion (1) . The stress tensor and heat-flux vector in this case contain turbulent con- 
tributions involving two- and three-point correlations, in addition to molecular 
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contributions. The turbulent correlations are unknown and must be approximated by a 
mathematical turbulence model. The turbulence models used in the code are called 
eddy-viscosity models and are represented in equations (6) by expressing y v and y e 
in terms of an eddy-viscosity function p-p, i.e.» 


y v = y + y T » 



( 8 ) 


where Pr.p is a turbulent Prandtl number. 

The dot product of the area vector and the flux tensor & appearing in 

equation (1) is the flux vector which can be expressed in terms of inviscid and vis- 
cous contributions as follows : 
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The viscous flux vector F^ 1 may be represented in matrix form by the sum 


= -We" 
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( 10 ) 


The above equations have been written explicitly for 2-D plane flow. For axi- 
symmetric flow, the corresponding equations can be obtained from the 2-D forms by the 
replacements 



(ID 


where H is a source function to be added to the right-hand side (RHS) of the v 
momentum equation, • 3 a u = S ax 3 a u S a y3 a v and r = y is the radius or distance 

from the axis of symmetry. 
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2.2 Finite Volume Method 

The finite-volume technique is used to express the equations of motion in dis- 
crete form. In this method, a lattice of mesh points is generated by specifying the 
Cartesian coordinates x and y as discrete functions of the computational coordi- 
nates K and n. In the discrete approximation, £ and n are given either integer or 
half integer values, and are used interchangeably with the subscripts i and j (see 
fig. 1). The discrete representation of the surface vectors and volume element is 
given by the following formulas : 

^nh.j-Ci/a) = -i(y i+i,j “ y i.j ) + " Xi >3 ) 

u i,j = (5 i+i " *i )(y j+i " y j> _ (x j+i _ 5 j )(y i+i ' y i } 

H “ J (x i,j + x i,j+i> ’ x j “ I (x i,j + x i+i.3 ) ’ etC 

In the axisymmetric case, the radial variable r *= y in equations (11) is rep- 
resented as an average of the four associated points in the expression for i), and an 
average of the two associated points in the expressions for t a . For the planar 
case, the surface vectors satisfy the numerical analog of the metric identity, equa- 
tion (4). A related identity is satisfied in the axisymmetric case. This feature 
leads to the important physical property of free-stream maintenance in which the dif- 
ference approximation of ^ vanishes when p,p,u are constant regardless of 

the difference approximation for 

In the finite— volume method the state vector is assumed to be an integrated 
average over the cell volume (area) . It is located at the centroid of a cell and is 
given integer subscripts, i.e., U. . . This notation is frequently contracted, e.g., 

1 9 J 

U., , = U. , . The surface vectors are located at cell boundaries and are given 

1+1, j 1+1 

mixed integer /half integer subscripts. 

2.3 Diagonal Transformation 

The Jacobian matrix of the inviscid flux vector F^ is defined by A a = 3F^/3U. 
This notation will be shortened here by dropping the subscript a. Thus, we write 
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pQ 


(PE + p)Q (13) 

mm mm 

Q = u • ? *» uS x + vSy 

In the following development, considerable use will be made of the diagonal rep- 
resentation of A. This is written 

A = SR _1 AR > 

A = diagCun^.Uj! + c,u^ - c) 

V ( 14 ) 

S = \t\ , n = ?/S , u n « u • ^ 

R = CNP , R" 1 « P" 1 N“ 1 C" 1 ^ 

The matrix R is a similarity matrix which diagonalizes A, SA is the diagonal 
matrix of eigenvalues of A, and is the component of velocity normal to the sur- 
faces defined by The matrix R is factored into three submatrices which are 


written 



( 15 ) 


where x = Y - 1 and the empty spaces in the matrices represent zeros. In the 
orthogonal matrix N, the elements and ny are the Cartesian components of 

i.e., n = Tn x + jny, and m = + jmy is a unit vector orthogonal to "£ with 

m^. - -n y , my = The component of velocity tangent to the surface ^ is given by 

the dot product u^j « u * Under the transformation R, the vector of conserved 
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variable differentials is transformed into a vector of differential Reimann variables, 
i.e. , dW = RdU or 
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(16) 


2.4 Time Differencing 

The basic algorithm used to solve equations (1) is an approximately factored ADI 
method based on the work of Beam and Warming (ref. 3). The implicit time-differenced 
equations are obtained by Taylor series expansions in t utilizing the relations 

3 t U = [U(t + At) - U(t)]/At = AU/At'' 

u(t + eat) ~ u(t) + eau 


F,J(t + eat) - F£(t) + 8A a (t)AU 
F£(t + 9At) - F£(t) + 6B a g(t)3 g AU J 


(17) 


In these expressions AU is the incremental or delta variable. At is the time step, 
and 6 is a constant (between 0.5 and 1.0). An approximate form of the expansion of 
is used in which the viscous matrices B a ^ are frozen in time. Substitution of 
the above expansions into equations (1) gives the unfactored form of the implicit 
algorithm 


[I + h'3 a (A a - B ag 3 g )]AU = -A3 a F a ’ 
h = At/u , h ? = 0h , a, 3 = 


(18) 


where the summation is taken on both a and 3> and I is the identity matrix. The 
value of 0 is normally taken to be 0=1 corresponding to first-order, backward- 
time differencing. 

The above equation is simplified by approximately factoring the implicit, or 
left-hand side (LHS) operators into £ and r\ operators. A further approximation 
consists in discarding the off-diagonal viscous matrices, i.e., B a g = 0, 3 £ a. By 
redefining B a = B a g, 3 = a, a chain of implicit factored equations in ADI form is 
obtained : 
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( 19 ) 


[I + h'a^(A^ - B^)]AU* = -h(3 5 F 5 + 3^) 
[I + h’8 n (A n - B n 3 11 )]AU - AU* 

U(t + At) = U(t) + AU 


By suitably defining the partial derivatives in terms of spatial difference 
operators, the above equations can be expressed as a system of coupled algebraic 
equations which may be solved by block— tridiagonal matrix inversions. This is a 
costly numerical process. A more efficient procedure, which is used in the code, is 
to simplify the left-hand sides of equations (19) utilizing the diagonal transforma- 
tion so that they become uncoupled and can be solved by scalar tridiagonal inversions. 
This process results in the diagonal algorithm (ref. 6). 

To develop the diagonal algorithm we write a typical member of equations (19) by 
dropping the subscripts on the matrices A and B, e.g., 

[I + h’3 ? (A - B3g) ] AU = AU C (20) 

The viscous matrix B is approximated by the identity matrix multiplied by the 
largest eigenvalue of B, i.e.. 


B = vl , 


v 


_S_ 

po 
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^v 



( 21 ) 


By introducing the diagonal decomposition of A, i.e., A = SR 1 AR, and making 
the further approximation of taking A, R, S, and v to be locally constant (so that 
they can be moved outside the partial derivatives) we obtain a diagonal, uncoupled 
system of equations or the diagonal algorithm 


[I + h"(A3 ? 


- vI3|)]RAU = RAU q 

h" « 6AtS/u 


} 


( 22 ) 


After making a suitable specification of the spatial derivatives in terms of finite 
differences, the above equations are solved recursively by scalar tridiagonal 
inversions . 

The approximations involved in the reduction to diagonal form diminish to some 
extent the time accuracy of the algorithm (ref. 6). For example, the nonconservative 
nature of the implicit operator causes some error in the propagation speeds of strong 
shock waves. However, the approximations do not affect the accuracy of a steady- 
state solution, which is governed by the spatial flux differencing of the right-hand 
side of equations (19) and which is normally taken to be second order. 
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2.5 Spatial Differencing 

To complete the specification of the numerical algorithm, the spatial deriva- 
tives in equations (19) and (22) must be replaced by finite differences. This is 
done by using central differences for the viscous terms and upwind differencing for 
the inviscid terms. 

Consider first the RHS explicit flux differencing of equations (19). To main- 
tain strict conservation regardless of the type of differencing used to approximate 
the fluxes, the difference approximation for ^F^ is written as 

a s F s = ( Vi+(i/2),j ~ ( Vi-(i/o,j (23) 

The (F ,-) / \ • are fluxes evaluated at the cell boundaries that must be expressed 
in terms of the state vector evaluated at the cell centers, e.g., U^j, 

l^i+i . , etc. To simplify the notation we will drop the subscripts t, and j from 
F^ in equation (23), i.e., 

= 3 £ F = F i+(i/ 2 ) ’ F i-(i/ 2 ) 

F i+(i/2) = ^i+(i/ 2 ) * *^i+(i/ 2 ) “ F i+(i/ 2 ) + F i+(i/ 2 ) 

(a) Inviscid flux differencing - The inviscid flux vector F i+(i/ 2 ) * s re P re ” 
sented in terms of a central average involving (U^) and = «^(U^ +1 ) and 

a dissipation function D^^^) : 



F i4-(i/2)„ “ 2 ^i+(i/2) * + t ^i+P “ l^i+(i/ 2 ) l D i+(i/ 2 ) ^ 

If Di+( 1 / 2 ) - 0, the differencing for 3^F reduces to pure second-order central 
differencing. 

The dissipation function is expressed in terms of special operators ^ and 6 
that are defined as follows: 

D i+d/2) - [R _1 (A^U +|A|^U)] i+(l/2) 

^>i+(l/2) = “« W i+(3/2) - SW i-(l/ 2 )> 

^)i+d/ 2 ) * B« W i+(x/2) - Y(fiM i+ ( 3 / 2 ) + W; 

(«W) i+( i/ 2) = iW i+(l / 2 ) = Ri+d/2) (U i+i “ U : 

The constants a, $, y control the accuracy of the spatial differencing for d^F’. 
Several cases are listed below. 



(26) 
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Spatial flux differencing options 


First-order upwind 
Second-order upwind 
Third -order upwind 
Second-order central 


3 “ 1 , a=Y=0 

3 = 1, a « y = 1/2 
3 = 1/3 , a = y = 1/6 
3 = free , a = 0 , y 


3/2 


(27) 


Most of the cases reported in this paper were done with the second-order upwind 
option. Assuming A is constant, spatial accuracy may be checked by Taylor series 
expansion of equation (24) using equations (25), (26), and (27). 

To completely specify the inviscid flux differencing, the method of averaging, 
i.e., computation of R i+( 1 / 2 ) in terms of u i an< * U i+1 , must be given. The type of 
averaging used here is the nonlinear averaging of Roe (refs. 7 and 8). The specific 
formulas used in computing ’&±+( 1 / 2 ) an< * ^i+(i/ 2 ) are 8*- ven below. 

a = (Pi+j/Pp 1 ^ 2 , H = E + p/p 

u i+(i/2) = (au i+i + u i )/(1 + a) 
v i+(i/a) " (av l+i + v i )/(1 + a) 

c l+(i/ 2 ) “ <r - lM(aH i+1 + H i )/(1 + a) - (u 2 +(l/2> + v 2 +(l / 2 ))/2) 



These formulas are used with equations (15) where it is noted that the unit-vector 
components n y , m^, being obtained from require no averaging. 

(b) Viscous flux differencing- As stated above, the viscous flux F M = ? • 8?" 
is centrally differenced using equation (24) . Explicit formulas for evaluating the 
heat-flux vector q at the cell boundary i + (1/2) are described below. 


^i+(i/2) “[u ( ^5 3 C e + s n 3 Ti e) ]i+(i/2) 


(29) 


The viscosity and volume element in this expression are evaluated by central 
averaging, i.e., 

f U h +(l /2) ' I <»i + »i«) (30) 

The surface vector (^)^ + ^ 1 y 2 ^ is given directly by the finite-volume construction. 
The 


normal derivative term is evaluated by 

(3 C e) i+(i/2) = e i+i " e i 


(31) 
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The two remaining terms are cross -derivative expressions in which a special averaging 
procedure is used to center known quantities at the cell face. These cross-derivative 
terms are given below. 

^i+(i/2) = 2 [( ^n ) j-(i/2) + ( ^p ) j + (i/ 2 ) ] i+(i/ 2 ) 

^n e ^i+(i/ 2 ) " 2" ^ e j+i "" e j-i^i+(i/ 2 ) 

The additional elements of , i.e., T and x • xj are determined by a similar 
construction. 

An option to use a simplified form of the Navier-Stokes equations called the 
thin-layer Navier-Stokes equations is included in the code. Under this approximation 
the derivatives tangential to a thin shear layer (say 3^) are discarded resulting in 
a simplified form of the equations. Use of this option gives a reduction of about 
15% of the computing time required for the full Navier-Stokes equations. 

(a) Implicit differencing- The differencing of the LHS or implicit terms of 
equations (22) follows a pattern similar to that used for the RHS or explicit terms. 
That is, upwind differencing is used for inviscid first derivatives and central dif- 
ferencing is used for viscous second derivatives. These differencing expressions are 
given as follows: 



where X is one of the eigenvalues u n ,u n ,u n + c,u n - c, and f is a corresponding 
element of the vector of differential Reimann variables, i.e., AW - RAU (see 
equation (16)). 

2.6 Boundary Conditions 

Boundary conditions must be specified for the flux vector F a and delta variable 
AW — RAU at the physical boundaries of the mesh. Unfortunately, space does not per- 
mit a detailed exposition of boundary techniques here; thus, only a brief summary will 
be given. A more detailed description is given in reference 1. 

We will discuss first boundary conditions for the inviscid Euler equations. 
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(a) Inviscid boundary conditions - The basic principle in specifying boundary 
conditions for the inviscid equations comes from the method of characteristics (MOC) . 
This principle states that the number of conditions or variables that may be arbi- 
trarily prescribed at a boundary is equal to the number of eigenvalues of A whose 
signs indicate signal propagation from the boundary to the interior domain. The 
remaining conditions or variables must be determined by some form of numerical extrap- 
olation from the interior to the boundary. This can be done in many ways. The 
method used in the code uses a set of compatibility relations associated with each 
eigenvalue and derived from a locally one-dimensional approximation. These compati- 


bility relations are listed below. 

Eigenvalue 

Compatibility relations 


“n 

V ■ 3 5 % “ 



“n + c 

3 £p + pc^Un - 0 , 

► 

(34) 

“n ” c 

3^p - pc3^u n » 0 J 




where s' = c/p^“ 1 )^Y i s an entropy variable and and u m are the normal and 
tangential components of velocity at the boundary. The compatibility relations are 
or are not used at boundary depending on the signs of the eigenvalues. The procedure 
is summarized below for a left boundary. 



Flow state 


Boundary 

r condition 

1 . 

Inflow, u n > 

0 

Prescribe 


2 . 

Outflow, u n < 

0 

Extrapolate 

3^ (s T »%) “ 0 

3. 

Supersonic, u n > 
inflow 

c 

Prescribe 


4. 

Supersonic, < 

outflow 

-c 

Extrapolate 

^(p,^) = 0 

5. 

Subsonic, [iijJ 

flow 

< c 

Prescribe 

Extrapolate 

f (p,^) = 0 
a^p - pca^Un = 


(a) Prescribe normal velocity 

(b) Prescribe pressure p 

(c) Prescribe total temperature c 2 + ^ ~T~~ ( u n + u m^ 
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The extrapolation conditions for 3^w are numerically approximated by either a 
first- or second-order form, i.e., 

3^w » 2(w ] . - w fi ) or 3w ] . - w I+i - 2w fi (35) 

where the indices I, I + 1, and B are indicated in figure 2. 

The above conditions result in the determination of a nonconservative state 
vector V = (s ' ,u m ,u n ,p) T at the boundary. The flux vector Fg, which is required 
at the boundary is constructed from V, i.e., Fg *= Fg(V). 

Inviscid boundary conditions on the delta variables AW — RAU are determined by 
a procedure similar to that described above. Details may be found in reference 1. 

Boundary conditions are needed on the differences ^j|.+(i/ 2 ) equations (26). 
The procedure presently used in the code for these terms is simple extrapolation, 
i.e., <SW I _( 1 / 2 ) = 6Wg - 6W i+ ( 1 / 2 ) (see fig. 2). This procedure may lead to some 
error, especially for inviscid flows. More rigorous conditions are described in 
reference 9 but these have not yet been incorporated into the code. 

(b) Viscous boundary conditions - For flows in which viscous effects are impor- 
tant at boundaries, the inviscid conditions described above must be modified. The 
basic principle in these cases is that the variables associated with the normal- 
velocity component (i.e., the tangential velocity u^ and the entropy s, or tem- 
perature T) may be prescribed independently of the sign of u^ at the boundary 
(since the equations are no longer hyperbolic but parabolic) . For viscous flows that 
are subsonic at the boundary, the pressure and normal velocity are determined by the 
same procedure used in the inviscid case. 

Additional boundary conditions may also be imposed. In the code, the option of 
using periodic boundary conditions (associated with airfoil "0" grids) is incorpo- 
rated. In addition, a centrifugal correction to the wall pressure (determined from 
the normal momentum equation) which is important for inviscid flows over highly 
curved surfaces, can also be used. 

2.7 Turbulence Models 

At the present state of development, there is no single model that holds a 
decisive advantage over all other turbulence models. For this reason several zero- 
and two-equation turbulence models are integrated into the code to investigate their 
impact on overall modeling accuracy and to facilitate further development. These 
models are described in detail in reference 5 and only a brief summary of them will 
be given here. The zero-equation model is a version of the Cebeci-Smith model 
(ref. 10) and is summarized below. 
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(a) Zero-equation model ( Cebeoi-Smith ) - 

y T * min(u TI ,y To ) , Pr T =0.9 
y T i 80 p£ 2 |^ x u| , y^ Q “ 0.0168pu e S* 

Z - 0.4yD , D = 1 - exp(“/p w T w y/26y w 

where y is the distance normal to the surface, u is the boundary-layer edge 

C 

velocity, <5* is the kinematic displacement thickness, and the subscript w in the 
expression for the damper, D, represents wall values. The above model is applicable 
to smooth-wall boundary layers, and a modified form must be used for trailing edge 
and wake flows (ref. 5). 

(b) Two-equation models - Two-equation models use two additional PDEs for vari- 
ables that are used to define the eddy-viscosity function, y T . An advantage of these 
models over zero-equation models is that they can be applied over a wide range of 
flow geometries without the need for a special treatment or tailoring of an algebraic 
length scale. The field equations for two-equation models can be expressed in the 
form 

o3 t ps + 8 a ? a . • (psu + 5 S ) - H g 

s ■ “ 1" 1 * ?s “ “ ir 
“T = S D !’ '‘s-*‘+5£ 




where and Pr s are modeling constants, k is the turbulent kinetic energy, 
co = e/k is the specific dissipation rate df turbulence, and e is the absolute dis- 
sipation rate. The vector H g represents a pair of turbulence source functions and 
D is a damping function the details of which are given in reference 5 along with 
and Pr s . With these models, a turbulent -pressure term, (2/3)pkI is normally 
included with the viscous stress tensor t of equation (6). 

The three two-equation models incorporated into the code are identified below. 


Two-equation models 


1. 

k - e 

Chien model (ref. 11): 

s 1 = k , 

S 2 = e 

2. 

k - m 2 

Wilcox-Rubesin (ref. 12): 

Si - k » 

s 2 = w 

3. 

q - u> 

Coakley model (ref. 5): 

s 1 «= viT , 

s 2 = U) 
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The upwind-flux differencing procedure described earlier for the Euler equations 
is extended in a straightforward manner to the convective flux of the model equations, 
i.e., G a = psu • $ a . The implicit factored equation for s (corresponding to equa- 
tions (22)) takes the (nonconservative) form 

[1 + bh M v + h”(Xa^ - v3|)]pAs - (pAs) Q (38) 

where X = and h" and v are given as before. The term involving the parameter 
b is a critically important implicit source term used to balance the explicit RHS 
source term H s . The importance and use of this term are discussed more fully in 
reference 5. 

The viscous no-slip boundary conditions used with the two-equation models are 
k * q - 0 and e = 0 for the k - e model, 3 w = 0 for the q - w model, and an 
analytic BC, w = 20]i/(0.15py 2 ) , for the k - u 2 model. At inviscid inflow and out- 
flow boundaries these variables are treated in the same manner as s and (see 
eqs. (34) and (35)). 

Optional boundary conditions using wall functions (ref. 13), or matching to the 
law-of-the-wall, are also provided in the code. Compared with the conventional pro- 
cedure, this procedure is applicable only for high-Reynolds-number flows (i.e., not 
transitional flows) but has the advantages of requiring fewer points to resolve the 
boundary layer and can lead to substantial reductions in computing time. Our experi- 
ence with this procedure is limited, however, and we will show no results of its use. 

2.8 Time Step Selection 

For time-accurate calculations, a spatially constant time step is used. For 
airfoils, the value of At chosen is At ~ L/200^, where L is the chord length of 
the airfoil and c^ is the free-stream sound speed. 

For achieving rapid convergence to a steady-state solution, a spatially varying 
time step is used. We define the local explicit time step in the a direction to be 

Ata or 

At a = u/(u • ? a + c|?J) (39) 

which is the characteristic time required for an inviscid disturbance to propagate 
across a cell at the velocity |u n | + c. Utilizing this characteristic time, two 
options for a spatially varying or local time step are incorporated in the code. 

They are 
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Local time-step options 


(a) At(£,n) - CFL • minCAt^.At^) 

CFL » O (5-20) 

(b) At(£,n) ■ CFL • max(At^,At^) 

CFL = O (1-10) 

In most of the computations reported in this paper, option (a) was used with a 
Courant number of CFL = 6-10. 

For viscous— flow problems at high Reynolds numbers where a finely spaced mesh is 
used adjacent to solid surfaces, the above formula (option (a)) leads to very small 
values of At in the neighborhood of the surface and results in slow convergence. 

In these cases, the formula is modified by placing a lower bound on At, i.e.. 

At > At mln - L/20c„ (40) 

This procedure results in a much faster convergence to steady state. 

In starting a calculation from a uniform initial state, a local time step using 
option (a) is normally used with CFL = 1. As the iterations progress, the CFL number 
is increased by a factor of 1.1 per step to its maximum value. For high— Reynolds 
number viscous calculations, the time-step limit of equation (39) is imposed gradually 
over a period of about 50 cycles . 

3. CODE ORGANIZATION AND DESCRIPTION 

3.1 Summary of TURF Subroutines 

The code is named TURF (turbulent upwind resolution of flow). A summaary of 
TURF subroutines is listed in table 1 along with a brief description of the function 
of each subroutine. 

Input data (in the form of data statements) must be specified in the routines, 
MAIN, MSHA, MSHB, INCA, OUTA, CPLOT, AND SPLOT. These data are extensive and will 
not be discussed here. A manual explaining the code and how to set up test cases is 
currently under preparation and will be available in the near future. 

3.2 Grid Generation 

Algebraic grid generation routines are self-contained within the code (i.e., 

MSHA and MSHB) . Examples of the types of grids that can be generated are shown in 
figure 3. They are basically of two types: (1) sheared parabolic and elliptic con- 

formal mappings and (2) generalized H-grids where the points on two arbitrarily 
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defined curves are joined by straight lines. Spacing functions using uniform, expo- 
nential, and cubic-spline spacing are used to space grid points in both coordinate 
directions (SPAC). Combinations of different gridding procedures can be used to 
build up composite grids such as the airfoil grid shown in figure 3(a). A routine 
to compute the coordinates of NACA four-digit airfoils is also provided. For other 
airfoils and shapes a general curve drawing routine (CURV) is used to generate 
numerical coordinates from a set of specific input coordinates using linear or cubic- 
spline interpolation. Two types of airfoil grids can be generated: C-grids and 

0-grids. To date most numerical experience has been with airfoil C-type grids. 

3.3 Data Loading System and Vectorization 

A special data loading and sweep system is used (in MAIN) that is designed to 
permit implicit integration across interior mesh boundaries (see fig. 4). This inte- 
gration is achieved through the use of two data-array systems (or commons) called BIG 
and SMALL. Each system contains the same variables (e.g., p, pu . . .) but while the 
array BIG contains the data for the whole grid, the array SMALL contains the data for 
a single ^ or ri grid line. Both array systems are one dimensional, as indicated 
below. 

(a) BIG array: RO(LV), ROU(LV), . . . 

LV « I + IDM * (J - 1), £ « I, n = J 

(b) SMALL array: R(L) , RU(L), . . . 

L - I or J 

where IDM is the maximum dimension in the I direction. The variables shown here 
are density RO and R, and u-momentum ROU and RU. 

The general procedure for performing difference operations (either explicit in 
RHS or implicit in LHS) along a particular £ or q grid line is to first load the 
corresponding data from the BIG array Into the SMALL array. The data in the SMALL 
array run sequentially while the data in the BIG array may not. Next, finite differ- 
ence operations are performed on the SMALL array variables. Finally data from the 
SMALL array are off-loaded back into the BIG array. In this manner a single subrou- 
tine is used for differencing in either direction, which simplifies and reduces the 
amount of coding required. This procedure also permits implicit integration across 
internal mesh boundaries (fig. 4), thus allowing larger time steps than explicit 
procedures . 

Except for the scalar tridiagonal solver routine (TRX), essentially all Impor- 
tant DO loops in the code are vectorized on the Ames CRAY XMP computer. Work is 
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currently under way to vectorize the TRX routine (by interchanging DO loops) and it 
is anticipated that when this is accomplished computing times can be reduced by a 
factor of 2 to 3. 

4 . RESULTS 

Representative results of inviscid and viscous flow calculations are shown in 
figures 5-9. These were steady-state calculations run with a spatially varying time 
step that used, except for one case, the second-order upwind method. Computing times 
on the Ames CRAY XMP computer for the viscous airfoil calculations using the full 
Navier-Stokes equations take 128 and 160 sec for the zero- and two-equation models, 
respectively, for 500 time steps using a 160 x 50 C-mesh. Steady-state lift and 
drag convergence is achieved between 200 and 400 steps. Computing times for airfoils 
using the Euler equations, on a 120 x 30 mesh take 48 sec for 500 steps and converge 
to a steady state in 100-300 steps. 

4.1 Inviscid Flows 

(a) Lifting transonic airfoil - Results of inviscid flow calculations are shown 
in figures 5 and 6. Computed surface pressures for the inviscid transonic flow about 
an NACA 0012 airfoil in free air are shown in figure 5. This figure compares results 
obtained using the present method with those obtained from reference 14 which used 
the central-differencing and flux-splitting (second-order upwind) methods. The 
three methods are in basic agreement except in the neighborhood of the shock wave 
where the present method shows a sharper shock capture. 

(b) Oblique shook reflection - Figure 6 shows a calculation of a shock-reflection 
problem investigated by Yee et al. (refs. 8 and 15). Figures 6(a) and 6(b) show con- 
tour plots and pressure distributions using the second-order upwind method. A uni- 
formly spaced 60 x 20 cell H-mesh was used in the calculations. Supersonic inflow 
and outflow boundary conditions were imposed at the upstream and downstream bound- 
aries, respectively, and free-slip (u^ = 0) conditions were given at the lower 
reflecting surface. The Rankine-Hugoniot oblique-shock relations were used to impose 
conditions along the top boundary of the mesh. 

Although the oblique shock wave is captured with reasonable accuracy, there is a 
distinct undershoot in the pressure distribution ahead of the shock wave and a 
smaller overshoot behind it. These may be essentially removed by a high resolution 
option (ref. 2) in which a switch to first-order spatial differencing is done at 
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locations of relative maxima and minima in the pressure distribution. Results using 
this option (labeled UW2II-HR) are shown in figures 6(c)^and 6(d). 

4 . 2 Viscous Flows 

(a) Lifting tvccnsonic aivfoit - Results of viscous turbulent flow calculations 
are shown in figures 7—9. Transonic calculations of an RAE 2822 airfoil are compared 
with experimental measurements (ref. 16) in figure 7. The calculations were done at 
the experimental geometric angle of attack and nominal Mach number and do not account 
for wind-tunnel-wall interference effects. 

The meshes used in the transonic airfoil calculations were 120 x 50 (coarse) and 
160 x 50 (fine) C-grids (see fig. 3(a)). The mesh points were exponentially spaced 
away from the airfoil surface and the cell spacing adjacent to the surface corre- 
sponded to a y + of 0.8. About 23 mesh points were contained in the boundary layer, 
and the outer boundary was placed at 10 chord lengths from the airfoil. The cell 
spacings in the streamwise direction over the central part of the airfoil were 
Ax/c = 0.036 and 0.024 for the coarse and fine meshes, respectively, with finer 
spacings used at the leading and trailing edges. 

Boundary conditions for the airfoil calculations were as follows: no-slip veloc- 

ity conditions at the airfoil surface, constant total enthalpy, entropy and tangential 
velocity component (u^) around the outer C part of the mesh, and constant (free- 
stream) static-pressure conditions along the vertical back (outflow) boundary. 

Computed Mach contours for the RAE airfoil, which illustrate the flow geometry 
and indicate overall computational quality, are shown in figure 7 (a) . Computed 
pressure distributions are compared with experimental measurements in figure 7(b). 

The computations were done with the q - to two-equation model and include results 
using both the coarse- and fine-mesh systems. Pressures obtained using the various 
other turbulence models showed very little difference from the pressures obtained 
with the, q - w model. Almost identical results were achieved using the coarse- and 
fine— mesh systems, except in the region of the shock wave where small differences can 
be observed. 

Skin-friction results are compared in figure 7 (c) . Since the computed results 
using the coarse mesh were very similar to those obtained with the fine mesh, only 
results obtained with the fine mesh are shown. For the zero-equation model, the 
boundary layer was assumed to be turbulent starting at the leading edge. The q - m 
two-equation model predicts (natural) transition at an x/c, 53 0.05. The other two 
models indicate fully turbulent boundary layers starting at the leading edge. Of the 
two types of models tested, the zero-equation model gives the best overall prediction 
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of skin friction downstream of the shock wave. This trend in model performance is 
usually reversed, however, for flows with extensive separation. An example of this 
is given in the oblique shock boundary-layer interaction flow to be discussed shortly. 

(b) Symmetric airfoil in a channel - Calculations of the transonic flow over an 
18%-thick circular arc airfoil in a channel are compared with measurements (ref. 17) 
in figure 8. The shadowgraph and experimental pressure distributions of figures 8(a) 
and 8(c) show the basic features of the flow — an oblique shock-wave pattern, a thick 
zone of separated flow, and nearly critical values of the pressure coefficient down- 
stream of the shock wave. 

The numerical grid used in these calculations was an 80 x 60 H— mesh, the middle 
segment of which is shown in figure 3(b). The top mesh boundary is fitted to the 
channel wall and a plane of symmetry is used so that only the flow in the upper half 
of the channel need be computed. Boundary conditions at the upstream boundary were 
taken as constant total enthalpy and entropy conditions. Constant pressure condi- 
tions were taken at the downstream boundary. Along the top wall of the channel 
inviscid free-slip conditions were used, and along the airfoil surface viscous no— slip 
conditions were used. Inviscid free-slip conditions were also used along the symme- 
try plane upstream and downstream of the airfoil. The turbulence models utilized in 

2 

the calculations included standard and modified forms of the Wilcox-Rubesin k - m 
model described in references 1 and 5. The numerical scheme for computing the invis- 
cid fluxes was an earlier version of the second-order central scheme reported in 
reference 1. 

Numerical calculations of the pressure distribution along the airfoil surface 
are compared with measurements in figure 8(c). A plot of computed Mach contours is 
shown in figure 8(b) which may be compared with the shadowgraph of figure 8(a). 

(c) Oblique shook boundary -layer interaction - As a final case, calculations of 
an oblique shock-wave boundary-layer interaction flow are compared with experimental 
measurements (ref. 18) in figure 9. The flow consists of an oblique shock wave 
impinging on a turbulent boundary layer at a free-stream Mach number of 2.9 and a 
unit Reynolds number of 5.7 x 10 7 . The shock wave is sufficiently strong to cause 
substantial separation with a corresponding plateau in the measured pressure distri- 
butions (fig. 9(c)). 

The numerical grid used in the calculations was a 100 x 50 H-mesh. Supersonic 
inflow and outflow conditions were used at the upstream and downstream boundaries, 
respectively, and a boundary-layer program was used to provide profiles of the state 
variables at the upstream boundary. Along the top surface, boundary conditions were 
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determined from the Rankine-Hugoniot conditions for oblique shock waves (in subrou- 
tine SHOCK) . Results using both the zero-equation Cebeci-Smith model and the two- 
equation q - a) model are shown. The second-order upwind scheme was used for the 
inviscid flux differencing. 

Computed Mach contours using the two models are shown in figures 9(a) and 9(b), 
and experimental pressure measurements are compared with computations in figure 9(c). 
It may be seen that the zero-equation model significantly underpredicts the amount of 
separation and the extent of upstream influence (indicated by the plateau in the 
experimental pressure distribution) . The use of the two— equation model substantially 
improves predictions of surface pressure. 

5 . CONCLUDING REMARKS 


In the preceding sections we have described a code for solving the compressible 
2-D and axisymmetric Navier-Stokes equations in arbitrary curvilinear coordinates 
using the finite volume approach. Important features of the code, the numerical dif- 
ferencing procedures, boundary conditions, turbulent models, and time-step selection, 
were discussed. In addition, the general layout and operation of the code, including 
self-contained grid generation and plotting routines, were described. Finally, 
applications to a representative set of inviscid and viscous flow problems were 
presented. 

Although much of the code is in a final state there are several areas in which 
improvements are needed. These are: (a) vectorization of the implicit operator, 

(b) time-step selection and convergence acceleration, (c) improved boundary condi- 
tions on the dissipation variables, <SW i+ ( x / 2 ), of equation (26), and (d) fully checked 
out law-of-the-wall boundary conditions for the turbulence models. Work in these and 
other areas is proceeding. 
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TABLE 1.- SUMMARY OF TURF SUBROUTINES 


Subroutine 

Function 

MAIN 

Input data, loading and time-stepping loops 

MS HA 

General-purpose grid generation 

MSHB 

Airfoil grid generation 

CURV 

Space curves 

SPAC 

Spacing functions 

FOIL 

NACA 4 digit coordinates 

MTRC 

Surface vectors and volume elements 

CURT 

Curvature of boundary surface 

MPLT 

Mesh plots 

INCA 

Initial conditions (uniform ICs) 

SHOCK 

State behind oblique shock wave 

PROF 

Initial starting profiles for boundary layer 

RHS 

Right-hand side control and flux differencing 

vise 

Molecular and turbulent viscosities 

FLX 

Flux vectors at interior cell boundaries 

RBC 

Flux vectors at external mesh boundaries 

SORC 

Turbulent source functions and time step 

WALL 

Law-of-the-wall boundary conditions 

DELT 

Boundary-layer thickness functions 

LHS 

Left-hand side control and update 

TRX 

Scalar tridiagonal solver 

LBC 

Implicit boundary conditions 

OUTA 

Output printer plots, field data, surface data 

MAPO 

Routine used for printer plots 

CPLOT 

Contour plots M, p, \p 

SPLOT 

Surface data plots Cp, Cf, <$*» 0 

HEAD 

Plot heading routine 
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Figure 1.- Finite volume representation. 
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Figure 2.- Boundary cell indexing. 
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Figure 3.- Representative mesh systems. (a) C-mesh system for lifting airfoils; 
(b) H-mesh system for circular-arc airfoil in a channel. 




Figure 4.- Augmentor wing mesh system and implicit integration path ( ). 
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Figure 5.- Comparison of computed surface pressures for inviscid transonic flow about 
an NACA 0012 airfoil: M = 0 . 75 , a » 2 0 , 80 x 20 cell C-mesh (79 x 31 mesh for 

ref. 14). 


METHOD UW2II METHOD UW2II HR 



U -o w 

x X 

PRESSURE DISTRIBUTIONS 

Figure 6.- Comparison of computed pressure contours and distributions for inviscid 
shock reflection problem: M = 2.9, 0 =29°, 60 x 20 cell H-mesh. (a) Contours, 

method UW2II; (b) pressure distributions, method UW2II; (c) contours, 
method UW2II-HR; (d) pressure distributions, method UW2II-HR. 
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MACH CONTOURS 



SURFACE PRESSURE 



SKIN FRICTION 



Figure 7.- Transonic flow about RAE 2822 airfoil, case 9: M « 0.73, Re = 6.5 x 10 6 , 

a = 3.19. (a) Computed Mach contours; (b) surface pressures; (c) skin friction, 

upper surface. 
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Figure 9.- Oblique shock-wave boundary-layer interaction: M= 2.9, Re/L « 5.7 x 10 7 . 

(a) Mach contours, zero -equation model; (b) Mach contours, q - to model; 

(c) surface pressures. 
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